//lc: size of mesh
lc = 4000;

//points
Point(1) = {-20000, -20000, 20000, lc};
Point(2) = {20000, -20000, 20000, lc};
Point(3) = {20000, 20000, 20000, lc};
Point(4) = {-20000, 20000, 20000, lc};
Point(5) = {-20000, -20000, 0, lc};
Point(6) = {20000, -20000, 0, lc};
Point(7) = {20000, 20000, 0, lc};
Point(8) = {-20000, 20000, 0, lc};
Point(9) = {-20000, -20000, -20000, lc};
Point(10) = {20000, -20000, -20000, lc};
Point(11) = {20000, 20000, -20000, lc};
Point(12) = {-20000, 20000, -20000, lc};

Point(13) = {-5000, 0, 0, lc/100};
Point(14) ={5000, 0, 0, lc/100};



//lines
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line(5) = {5,6};
Line(6) = {6,7};
Line(7) = {7,8};
Line(8) = {8,5};
Line(9) = {9,10};
Line(10) = {10,11};
Line(11) = {11,12};
Line(12) = {12,9};
Line(13) = {1,5};
Line(14) = {2,6};
Line(15) = {3,7};
Line(16) = {4,8};
Line(17) = {5,9};
Line(18) = {6,10};
Line(19) = {7,11};
Line(20) = {8,12};

Line(21) = {13,14};




//faces 围成面的线条有方向
Curve Loop(1) = {1, 2, 3, 4};
Curve Loop(2) = {5, 6, 7, 8};
Curve Loop(3) = {4, 13, -8, -16};
Curve Loop(4) = {14, 6, -15, -2};
Curve Loop(5) = {13, 5, -14, -1};
Curve Loop(6) = {16, -7, -15, 3};
Curve Loop(7) = {9, 10, 11, 12};
Curve Loop(8) = {17,-12, -20, 8};
Curve Loop(9) = {18, 10, -19, -6};
Curve Loop(10) = {17, 9, -18, -5};
Curve Loop(11) = {20, -11, -19, 7};


Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};
Plane Surface(5) = {5};
Plane Surface(6) = {6};
Plane Surface(7) = {7};
Plane Surface(8) = {8};
Plane Surface(9) = {9};
Plane Surface(10) = {10};
Plane Surface(11) = {11};


//Vol
Surface Loop(1) = {1,6,3,5,4,2};
Surface Loop(2) = {2,11,8,10,9,7};

Volume(1) = {1}; 
Volume(2) = {2};

//Add 3 physical tag
Physical Volume("Air",1) = {1};
Physical Volume("Ground",2) = {2};


Field[1] = Distance;
Field[1].CurvesList = {21};
Field[1].Sampling =100;

Field[2] = Threshold;
Field[2].InField = 1;
Field[2].SizeMin = lc/200;
Field[2].SizeMax = lc;
Field[2].DistMin = 100;
Field[2].DistMax = 10000;

Background Field = 2;

// 外边界各面（示例命名，可按需改名）
Physical Surface("Outer") = {1,3,4,5,6,7,8,9,10,11};
// 如需空气-地层界面那一张中间面，也可以：
//Physical Surface("Interface")  = {2};

//Field[1] = Box;
//Field[1].Thickness = 1000;
//Field[1].VIn = lc/100;
//Field[1].VOut = lc;
//Field[1].XMin = -2000;
//Field[1].XMax = 2000;
//Field[1].YMin = -2000;
//Field[1].YMax = 2000;
//Field[1].ZMin = -1500;
//Field[1].ZMax = -500;

//Background Field = 1;

//Field[1] = Distance;
//Field[1].CurvesList = {33};
//Field[1].Sampling =100;

//Field[2] = Threshold;
//Field[2].InField = 1;
//Field[2].SizeMin = lc/400;
//Field[2].SizeMax = lc;
//Field[2].DistMin = 50;
//Field[2].DistMax = 20000;

//Field[3] = Box;
//Field[3].Thickness =200;
//Field[3].VIn = lc/100;
//Field[3].VOut = lc;
//Field[3].XMin = -2000;
//Field[3].XMax = 2000;
//Field[3].YMin = -2000;
//Field[3].YMax = 2000;
//Field[3].ZMin = -1500;
//Field[3].ZMax = -500;

//Field[4]  = Min;
//Field[6].FieldsList = {2,3};

//Background Field = 4;

//Field[2] = Box;
//Field[2].Thickness = 100;
//Field[2].VIn = lc/400;
//Field[2].VOut = lc;
//Field[2].XMin = -500;
//Field[2].XMax = 500;
//Field[2].YMin = -500;
//Field[2].YMax = 500;
//Field[2].ZMin = -50;
//Field[2].ZMax = 50;


//Mesh.MeshSizeExtendFromBoundary = 0;
//Mesh.MeshSizeFromPoints = 0;
//Mesh.MeshSizeFromCurvature = 0;

// This will prevent over-refinement due to small mesh sizes on the boundary.

// Finally, while the default "Frontal-Delaunay" 2D meshing algorithm
// (Mesh.Algorithm = 6) usually leads to the highest quality meshes, the
// "Delaunay" algorithm (Mesh.Algorithm = 5) will handle complex mesh size
// fields better - in particular size fields with large element size gradients:

Mesh.Algorithm = 5;

